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THE RELATION OF FINITE ELEMENT AND FINITE DIFFERENCE METHODS 

By 

Marcel Vinokur 
SUMMARY 

Finite element and finite difference methods are examined in order to 
bring out their relationship. It is shown that both methods use two types 
of discrete representations of continuous functions. They differ in that 
finite difference methods emphasize the discretization of independent 
variables, while finite element methods emphasize the discretization of 
dependent variables (referred to as functional approximations). An 
important point is that finite element methods use global piecewise 
functional approximations, while finite difference methods normally use 
local functional approximations. A general conclusion is that finite 
element methods are best designed to handle complex boundaries, while 
finite difference methods are superior for complex equations. It is also 
shown that finite volume difference methods possess many of the advantages 
attributed to finite element methods. 


INTRODUCTION 


The theoretical prediction of a three-dimensional flow past an arbitrary 
body requires a numerical solution. The traditional approach, which has been 
highly developed, is to use a finite difference method. Recently, the finite 
element method has been proposed as an alternative procedure. In order to 
evaluate the relative advantages and disadvantages of the two methods, it is 
essential to understand their common basis as well as their fundamental differ- 
ences. The present work attempts to do this by showing that both methods use 
two types of discrete representatations of continuous functions. The differ- 
ences in the two methods stem from the relative emphasis given to these 
representations . 

Since both finite difference and finite element descriptions employ 
different notations, each with a myriad of indices, we will use a rather 
cavalier notation with a minimum of indices. The notation, as well as some 
mathematical concepts that may be somewhat unfamiliar, are discussed in the 
next section. This is followed by a description of the two types of discrete 
representations of continuous functions. This framework is then used to 
examine and relate the finite difference and finite element methods as applied 
to continuous field problems. 

MATHEMATICAL PRELIMINARIES AND NOTATION 
A lower case letter will denote a function of real variables, e.g., 

u = ffx) . (1) 

The letter f denotes the functional rule, while x stands for a set of independ- 
ent variables (which may be general, curvilinear coordinates) spanning a 
domain V with a boundary S. The dependent variable u can represent a vector 
set of unknowns, in which case (1) is a set of equations. If the dependence 
on only some of the independent variables will be discretized, we will write 
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(ID as 


u = f(x,t) , (2) 

where the dependence on the variables x wili be discretized, while the 
dependence on the remaining variables t will remain continuous. A subscript 
will denote a partial derivative, or a component of a gradient. Thus, the 
differential of (1) is written as 

du = u x dx = f x dx , (3) 

where a summation or dot product is implied. On the other hand, a divergence 
will be denoted by the symbol 3/3x. Integrations over a domain or a boundary 
will be indicated by the letters V or S under the integral sign. If n is 
the normal at the boundary, the divergence theorem can be written as 

.' || d» - | fnto . (4) 

V s 

(Note that the sumbol dx has three different meanings in (3) and (4)). 

A capital letter will denote an operator acting on a set of functions. 


e.g. , 


v = F(u) . . (5) 

Here F is the operator rule, u(x) stands for a set of functions, and v(x) the 
resulting function(s) after performing the operation(s) . A local operator 
involves only algebraic or differential operations, while a non-local operator 
involves shifting or integral operations. Let 6u(x) denote the variation of 
the function u(x), which is a small change in u(x), keeping x fixed. Thus, 


6x = 0 , (6) 

and from their definitions, variation and differentiation are commutative, i.e.. 


Su x = (6u) x . 

The "differential 11 of (5) is then defined as 


(7) 
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6v = <$F(u, u) = lim F(u+6u) - F(u) , (8) 

IIM l * ° 

where | |<5u| | is some norm measuring the magnitude of 6u. Here 6F(u,<5u) is 
called the Frechet differential of F at the point u in the direction 6u. If 
F is a local operator, one can define a "derivative". Consider the operator 

F (u) = f(x,u,u x ,u xx ) . (9) 

Using (6) , (7) , and (8) , one can write 

<5F(u,6u) = f u 6u + f u (<Su) x + f u (6u) xx . (10) 

Using the rule for product differentiation, (10) can be transformed into 


6F(Mu) = {f u - & [f u - & (f u )]W» ♦ 4 { t f u - -k (f u ^ 6u + f u 6u x } ‘ 


XX 


XX 


By analogy with (3), the operator 


F u = f u ■ [f u x " 37 


( 11 ) 


( 12 ) 


can be called the Frechet derivative of (9) at the point u. In general, the 
the Frechet differential of any local operator F(u) can be expressed as 

6F(u,6u) = F u 6u + (F q 6u + Fjfit^ + F 2 6u xx + ...) . (13) 

In order to determine if a given operator is a Frechet derivative, one must 
define an adjoint operator. If 6u 1 and <$u 2 are two arbitrary variations of u, 
one can obtain from (9) the expression 


6u 2 5F(u,6u 1 ) = {f u <5u 2 - 3- [{f u - (f u )}6 u 2 - f u «u 2x ]}6u 1 

X XX XX 

* h { ' f u - K «„ * f „ 5 “ 2S V - f „ 5u 2x 5 V ' 

X XX XX XX 


(14) 


The coefficient of 6u^ in (14) has the form of a Frechet differential of some 
other operator in the direction <Su 2 . We can thus define an operator F(u) 
(within an arbitrary additive function of x) which is adjoint to F(u), such 
that 
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( 15 ) 


. f u Su - 4 U - ^ (f u )]8» - f tu >. 

X XX XX 

In general, for any local operator F(u) we can write 

6u 2 6F(u,6u 1 ) = 6u 1 6F(u,6u 2 ) + ^-(Foo 6 ^ 611 ! + F 01 6u 2 5u lx + F 10 (Su l <Su 2x + '’ ’ *^ 16 ^ 

An operator is self-adjoint if F (u) = F (u ) . Given the operators G(u), Gq(u), 
G-^fu), etc., the conditions under which 

G6u + (Gq6u + G^6u x + . ..) 

is equal to a Frechet differential 6F(u, 6u) can be easily determined from 
the relation 

dFfujdu^^ + 6 u 2 ) = dFfu^u^) + 6F(u + 6 Uj,6u 2 ) = 6F(u,6u 2 ) + 6F(u 1 + 6u 2 ,6u^). 

(17) 

The condition is found to be 

fiu^Gfujdup + ^ [dG^u^upd^ + 6G 1 (u,6u ]L )6u 2x + ...] 

3 (1®) 
= du^Gfu ^u 2 ) + ^ [6G 0 (u,6u 2 )6u 1 + 6G 1 (u,6u 2 )6u lx + ...]. 

It follows that G(u) must be a self-adjoint operator. If (18) is satisfied, 

one can easily show that F(u) is given (within an arbitrary additive function 

of x) by ^ ^ ^ 

F (u) = u j G(Xu)dX + [ f G 0 (Xu)dX + Ujc f G^XiOdX +...]. (19) 

0 0 0 

An operator acting on a set of functions which results in a real number is 

called a functional. (The norm | | 6u | | in (8) is an example.) The discussion 

of Frechet differentials and adjoint operators reveals the presence of annoying 

divergence terms. Since these can, in a sense, be removed using the divergence 

theorem (4), this suggests that a useful functional is the integral of an 

operator over the domain of x, i.e., 

I(u) = J F(u)dx. C20) 

V 

Using (13) and (4), it follows that 
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6l(u>6u) = | 6F(u,6u)dx = J F u 6udx + J (nF Q 6u + nF^i^ + nF 2 5u xx + , ..)dx . 


( 21 ] 


Expressing the gradient at the boundary in terms of normal derivatives, this 
can be written as 


61 (u,6u) = J F u 6udx + J (F Q 6u + F^i^ + F^u^ 


+ . ..)dx , 


( 22 ) 


where the subscript n signifies a normal derivative, and F Q , F^, ¥^, etc., are 
again newly defined operators. Equation (22) is the basis for a variational 
principle, which is the starting point for one form of the finite element method. 


DISCRETE REPRESENTATIONS OF CONTINUOUS FUNCTIONS 
Given an arbitrary function of the form (2), the most direct way to 
discretize the dependence on the variables x is to discretize x itself. The 
simplest procedure is to choose a set of N arbitrary points x^(i = 1 to N) , and 
to specify an approximation to u at those points. We thus define N. functions 
of t, 

u?(t)^f(x i ,t) , (23) 

where the superscript * signifies an approximate representation. We will refer 
to this as a Lagrange representation. In finite element terminology the points 
x^ are called nodes, and the functions u^(t) are sometimes called nodal 
parameters. A more sophisticated procedure, requiring a smaller number of 
points, is to specify also approximations to derivatives of u (which in the 
most general case need not be consecutive). An example would be to specify 
the set of first partial derivatives, 

u xi (t ) ^ f x ( x i^) • t 24 ) 


Such a representation will be called Hermite. Note that each point (node) 
would now have associated with it more than one parameter. If u represents 
a set of dependent variables, it is possible to discretize each by a different 
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set of discrete points. This is often done in practice. 

An alternative procedure is to divide the domain of x into N discrete 
volume elements V 1 (i = 1 to N) , and to specify an approximation to some 
functional of u defined over V 1 . A typical choice would be the integrated 
average 

u 1 (t) f(x,t)dx . (25) 

V 1 

By analogy with a Hermit e representation for point discretization, we can 
define higher order representations for volume discretization by specifying 
approximations to integrated higher moments of u. Volume discretization is 
useful in the finite difference solution of equations written in divergence 
(conservation) form. It is also necessary to define piecewise functional 
approximations (see below). In finite element terminology, the volume elements 
V 1 are called finite elements. 

Both types of discrete representations involve two degrees of freedom. 

One is the arbitrariness in the location of the points x^ (or volume elements 
V 1 ) . Any knowledge about the behavior of the function to be approximated can 
be used to make a judicious choice. The other freedom is the choice of the 
number and nature of parameters to specify at each point (or volume element). 
Here the nature of the equations and the numerical scheme can be a determining 
factor. 

Discretization of Dependent Variables 

The point discretization discussed above cannot represent integrals, or 
derivatives of higher order than the order of the representation. Also, a 
given representation gives no direct information at points other than the 
discretization points. Therefore, in order to obtain a numerical solution, 
one must also utilize (even if implicitly) an analytic representation of the 
arbitrary function. Any analytic function can be represented as an infinite 
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series in a complete set of chosen functions (providing the series converges). 

An obvious discretization is to choose N terms, and let the coefficients be 
the discretization parameters. We will generalize this notion, and approximate 
u by 

u (x,t) ** g[x; c i (t)] , (26) 

where g is any arbitrary, chosen function of x and the N parameters c^(i = 1 to 
N) . The parameters c^ are themselves functions of the undiscretized variables 
t. If u stands for a set of dependent variables, each one can be represented 
by a different function g, and the parameters c^(t) would be sets of 
parameters . 

A general representation which is nonlinear in the c^ cannot be easily 
integrated, and differentiation can rapidly lead to very complex expressions. 

For this reason, it is normally used only in curve fitting, and to approximate 
purely algebraic terms. An exception is the rational function approximation 

* C 00 (t)+ C 01 (t)x + + ••• 

u(x,t)**-2° — — j * (27) 

c 10 (t) + c n (t)x + c 12 (t)x + ... 

whose derivative maintains a simple form. Since (27) has some advantages over 

a polynomial, it has found uses in solving equations involving only local 

operators. In general, though, one chooses a linear representation in the c^, 
of the form 

* N 

u (x,t) ** T. c. (x) , (28) 

i=l 1 1 

where the <f^(x) are an arbitarily chosen set of linearly independent functions, 
sometimes referred to as basis functions. Since the basis functions should 
be easily integrated and differentiated, they are often taken to be powers of 
x, so that (28) becomes a polynomial in x. Other popular choices are 
trigonometric and exponential functions. Representations (26) and (28) will 
be referred to as functional approximations, or approximation by trial functions. 
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An important specialization of the linear representation (28) is to 

★ 

combine it with the point discretization (23) by requiring that u equals the 
* 

nodal parameters u^ at a set of N nodes x ^ , i.e.. 


u. (t) c.(t)<f>. (x.) . 

3 i=l 1 1 3 


(29) 


Since the <f>^(x) are linearly independent, one can always choose a set of x^ for 

which the matrix ^(x^) is non-singular, and thus solve for the ^(t) in terms 

* 

of the nodal parameters u^.(t). The c^ (t) are then said to be determined by 
interpolatory constraints, and the approximation (28) is then called a 
Lagrange interpolate to f(x,t) at the nodes x ^ . It can be represented 
directly in terms of the u^ (t) by introducing new basis functions <f>^(x), 
called canonical basis functions, with the defining property 


f>. (x.) = 6. . , 
i 3 ij 


(30) 


where 6^ is the Kronecker delta. They can be easily obtained from the original 

basis functions <f>^(x) by seeking the representation 

_ N 


k(x) = E c..ef> .(x) . 
3 i=l 31 1 


(31) 


It follows from (30) and (31) that 


N 

Z 

i=l 


c jiW 



(32) 


Since ^(x^) is non-singular, the coefficients c^ are uniquely determined by 
(32). The Lagrange interpolate to f(x,t) at the nodes x^ can thus be expressed 
succinctly as 


* N * ^ 

u (x,t) ^ L u.(t)£(x) . (33) 

i=l 1 1 


Canonical basis functions can also be defined for Hermite interpolation. 
For a first order representation, defined by (23) and (24), one can introduce 
the functions (x) and(f> il (x), satisfying 
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(34) 


^iO Cx j 5 = 6 ij * ^iox^j 5 = 0 

and 

^(Xj) = o , ^ ilx (x.) = «... . (35) 

These can be obtained in a manner analogous to that described above for Lagrange 
canonical basis functions. The Hermite inteipolate to f(x,t) at the nodes x i 
can then be written as 

u (x,t) [uT(t)£. 0 (x) + u^ftD^jCx)] , (36) 

i 

where the summation is over the total number of nodes. More general Hermite 
interpolates can be similarly formed. 

Piecewise Functional Approximation 

A single representation of the form (26) or (28) will be poor approxi- 
mation for functions that undergo rapid variation in the x domain. It is 
also difficult to construct such representations for domains with complex 
boundaries when the x domain is multidimensional. It is then advantageous 
to combine such representations with a volume discretization, and define a 
separate representation, in each of M volume elements V* 1 , of the form 

u (x,t) ^g^[x;c|(t)] (xeV^, i=i to N^) (37) 

in the general case, or 

*- N-* . 

u 3 (x,t)^ Z ckt)^(x) (xeV 3 ) (38) 

i=l 1 1 

in the linear case, where N- 1 is the number of parameters in element V- 1 . Such 

a representation is called a piecewise functional approximation, or approxi- 

*i 

mation by piecewise trial functions. If the approximations u J (x,t) are 
independently chosen in each volume element, the resulting global representa- 
tion would be discontinuous. 
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A representation with some degree of continuity requires matching 
conditions at interelement boundaries, which effectively limits one to the 
linear case (38). A practical method is to determine the c?(t) by interpo- 
lator constraints. We thus superimpose on the volume discretization an 
independent point discretization defining a set of N nodes x^ and associated 
nodal parameters. Matching is simply obtained by locating some of the nodes 
on interelement boundaries, where they are shared by more than one element. 

The N* 1 nodes belonging to element therefore satisfy the inequality 

M 

Z N J > N . (39) 

j=l 

One can again choose the set of nodes x^ so that ^(x^) is nonsingular 
(x^e V^) in each element . This condition will be sufficient to obtain 
continuity for an arbitary set of <f>^(x) if x is a one- dimensional variable, 
but continuity for multidimensional domains imposes restrictions on the set 
<{^(x). To show this clearly, we first discuss the one-dimensional case, but 
in a manner that can be immediately generalized to several dimensions. 

One-dimensional Representation . Let x be one -dimensional, and consider a 
piecewise representation (38) that is everywhere continuous, but whose 
derivatives can be discontinuous fct interelement boundaries. It is therefore 
sufficient to choose Lagrange interpolation, placing one node at each inter- 
element boundary, and additional nodes in the interior of each V*^ for which 
N-* > 2. One can then again introduce new basis functions (x), called 
Lagrange cardinal basis functions, satisfying 

- % < 4 °> 
for all i and k. In those elements V J which do not contain x^ 4 ) i (x) must 
interpolate to zero at all the element nodes. Since <f>?(x k ) is assumed non- 
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singular, it follows that ^(x) = 0 in those elements. Thus ^(x) is non- 
zero only over those elements containing node x i , i.e., two adjoining elements 
for a boundary node and a single element for an interior node. One-dimensional 
Lagrange cardinal basis functions for elements containing one interior node 
are sketched in the top row of figure 1. Since the interelement boundaries 
consist of one point at which an interpolating node is located, the functions 
*• (x) are continuous. Consequently, the global representation 

* N * r- 

u (x,t) ^ Z u (t) *.(x) (41) 

i=l 1 x 

is also continuous. The localized nature of the <P i (x) has important computa- 

* 

tional advantages. For example, integrals of products of u (x,t) over the 
domain define matrix elements k^ given by 

k ij = . (42) 

V 


It is seen that k.^ = 0 unless nodes i and j are contained in the same 

element, so that k^ is a sparse matrix. Similar results hold for integrals 

* 

of products of derivatives of u (x,t). 

In finite element applications it is convenient to define for each 
element a set of element cardinal basis functions ?^(x) satisfying 

1^(x k ) = 6 ik (Xj^evJ) . (43) 


If one extends the^(x) by defining them. to equal zero if x^ or x lie outside 
of v\ i.e.. 


^(x) =0 if x^^ or xeV^ , 

* j 

one can then represent u (x,t) for each V J as 

u **(x,t) ** E u. (t)$j (x) = H . u i (t)^ (x) , 
i=l ieV J 

where the second form results from (44) , It also follows from (44) that 


(44) 


(45) 
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(46) 


u i (x,t) =0 if x£V** . 


Using (43) through (46), one shows immediately that global and element 
representations are related by 



11 

<J>. (x) = Z & (x) 

3-1 

* M *• 
u (x,t) = 2 u 3 (x,t) . 

(47) 

and 

(48) 


3=1 



In order for (47) and (48) to be valid at interelement boundaries, the volume 
elements V* 1 must be considered disjoint, and to butt together at the 
boundaries. Element basis functions <H(x) are sketched in the bottom row 
of figure l. 

Another useful computational device is to define for each element a 

set of local coordinates x^ , each related to the global coordinates x through 

transformations x = x(x^) and x J = (x) . (A special case is x^ = x) . 

The nodes contained in each element V** can then be designated as x^ , where 

i is a local node number (i = 1 to N^), completely independent of its global 

node number. Thus there exist mapping relations which map local node numbers 

into global node numbers, and vice versa. The local nodal parameters 

attached to a local node x^ are designated as u^ J (t). The element cardinal 

basis function corresponding to local node x? would then be written as 

^(x j ), and the representation of u (x,t) in V** becomes 

. . 

u J (x 3 ,t) = £ u . 3 (ttf? (x a ) . (49) 

i=l 1 1 

The use of local coordinates can result in functions ^(x* 1 ) that are easy to 
manipulate analytically. A major advantage results if all the volume 
elements are geometrically similar in x space (which is trivially so in 
one dimension), since then they can all be described by the identical local 
coordinates. If the same set of basis functions ^(x*) is chosen for each 
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element, and the nodes x? are defined at geometrically similar locations, 
the element cardinal basis functions ??(x J ) will also be the same for all 
elements. It is thus possible to create a single subroutine, valid for all 
elements, in order to perform calculations for a single element. Of course, 
in summing the results to obtain a global solution, the coordinate trans- 
formations and node number mappings must be invoked. 

If continuity of derivatives is required for the piecewise representa- 
tion, one must use Hermite interpolation. It is only necessary to define 
derivatives at boundary nodes, and not at interior nodes. In fact, in 
most applications of piecewise Hermite interpolation, nodes are only 
defined at interelement boundaries. The extension of this subsection to 
piecewise Hermite interpolation follows the general manner indicated by (34) 
through (36) for the case of a single, global Hermite interpolation. 

Tensor Products 

A piecewise representation can be easily obtained in several dimensions 
if the global boundaries of the domain lie along the coordinate surfaces. 

One can then choose volume elements and nodes to lie along coordinate 
surfaces, and construct cardinal basis functions which are products of one- 
dimensional cardinal basis functions known as tensor products. We indicate 
the process for two dimensions, departing from our notational convention, 
by using x and y to represent the two (not necessarily Cartesian) 

coordinates. 

k ^ 

Let V , x^, and <J)^(x) be one-dimensional volume elements, global nodes, 
and Lagrange cardinal basis functions along the x coordinate. Similarly, 

% ft* 

define V , y^ , and <Pj (y) to be one-dimensional volume elements, global nodes, 

and Lagrange cardinal basis functions along the y coordinate. These define 

k& 

two-dimensional volume elements designated as V , and the double index node 
number ij for the node located at and y .. The function 
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* i; j(x,y) = <f>. (x)4>j (y) 


(50) 


has the property 


^ij ** x m ,y iP ~ ^im^jn 


(51) 


and is therefore a two-dimensional Lagrange cardinal basis function. 


Consequently, 


where 


N N 

3d y ^ 

u*(x,y,t) Z l u* . (t)f. . (x,y) 
i=l j-1 13 13 


(52) 


u ij (t) ^ fCXi./j.t:) (i=l to N x , j=l to Ny) . (53) 

Representation (52) is everywhere continuous, and can be extended to higher 
dimensions and to the case of Hermite interpolation. 


General Multidimensional Representation 

If the global boundaries of a multidimensional domain are too complex to 
allow a tensor product piecewise representation, one must use volume elements 
of a more general shape. All the results of the subsection on the one- 
dimensional representation can be immediately generalized, with the exception 
of the continuity conditions. If x^ is an interior node located in element 
V* 1 , we require that ^(x) (or *$?(x)) equals zero on the boundaries of V** . 

But this is only guaranteed at the boundary nodes of V* 1 . Thus the 
combination of basis functions <f^(x) in (38) and nodes x^ cannot be 
arbitrarily chosen, but must be such as to yield ^(x) = 0 on the boundary 
for interior nodes of . If x i lies on one or more boundaries of V* 1 , then 
we require that??(x) = ?^(x) on each boundary for all other volume elements 
sharing that boundary. In addition we still require that cj^(x) = 0 on 
the boundaries of that do not contain x^ Piecewise Hermite interpolation 
puts even more stringent requirements on the <f>?(x). 

Up to this time, the shape of the volume elements V** and the nature of 
the basis functions <f>?(x) have been considered arbitrary. The above-mentioned 
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continuity requirements effectively limit one to triangles (or tetrahedrons) 
in x space (which could be curvilinear in physical space), and polynomials 
in x for the <f>?(x). It also puts restrictions on the location of the nodes 
x i « The simplest case is Lagrange interpolation with linear basis functions 
<t>?(x), for which one only requires nodes at the vertices of the triangles 
(or tetrahedrons) . The finite element literature is replete with various 
combinations of nodes x^ and corresponding cardinal basis functions 
(usually called shape functions) ^(x), for both triangles and tetrahedrons, 
and for Lagrange and Hermite interpolation. 

The polynominal nature of the <f>-?(x) also allows one to estimate the 
* 

errors in u (x) , when f(x) is assumed exact at the nodes (so that (23) is 
an exact equality). (We suppress the dependence on t for the moment.) Such 
estimates are derived in reference 1, where it is shown that the error bound 
for Lagrange interpolation over a triangle is inversely proportional to the 
sine of the smallest angle. This would rule out extremely acute triangles. 
Actually, the author has shown (ref. 2) that the sine of the largest angle 
enters into the error bound, ruling out only extremely obtuse triangles. 

For the simple linear case, the author obtained least upper bounds for the 
errors. Let 


M E l f xxl max 

be the maximum absolute value of the second directional derivative of f in 
any direction, at any point in the triangle. If 6 and h are. the maximum 
angle and side of the triangle, then the results are 


2 

— 


and 


|u*-f| <“■ 

I 1 sir • 


where | fu -£) | is the magnitude of the gradient of u -f . 


(54) 


(55) 

(56) 
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Since arbitrary, curved global boundaries cannot be easily fit by a 
global curvilinear coordinate, one is faced with the need to use curvilinear 
elements in x space. This can be done if one can find transformations x(5) 
which transform the curvilinear elements in x space into straight sided 
"parent" elements in £ space. The representation of u over the curvilinear 
element is only approximate, being accurate only at the nodes. It is there- 
fore sufficient to treat the transformations x(£) in the same manner. This is 
the basis for the isoparametric transformations developed by Irons (ref. 3). 

Let x^ = x be the local coordinates for the curvilinear element, and x? be a 
set of local nodes chosen on the boundary of (and possibly within) the 
element. (There must be at least 3 nodes per side and at least 4 nodes per 
face to define curved boundaries.) In the transformed 5 plane, let £? ,5 ^ , 
and *$(£**) be local coordinates, local nodes, and element cardinal basis 
functions for the corresponding "parent" element. Then the isoparametric 

transformation has the approximate representation, 

* . . . . . 

x J (5 J ) ^ 2 x 3 VAZ 3 ) . (57) 

i=l 

(In some cases it is practical to use a lower (higher) number of nodes and order 
of basis function to represent the geometric transformation than are used to 
represent u J (x) over the element. Such transformations are then called sub 
(super) parametric.) Using (57) and its derivatives, integrals over element 
in x space can be transformed into integrals over the "parent" element in 5 
space. 

Splines 

The piecewise representations discussed so far involved only interpolatory 
constraints to determine the c?’(t) in (38). If additional smoothness constaints 
are imposed by matching higher derivatives (than those prescribed by interpola- 
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tion) at boundary nodes, the representations are called splines. The additional 
continuity requirements on the 4>?(x) are so great for general multidimensional 
elements as to make such spline representations totally impractical. We are 
thus restricted to one-dimensional splines (and their tensor product generaliza- 
tion to higher dimensions). In practice, splines are further limited to 
volume elements with Lagrange interpolatory nodes only at the two ends of each 
element. Thus for a division of the one- dimensional x space into M volume 
elements, the total number of nodes N = M + 1. Smoothness constraints are 
applied at the M - 1 nodes that lie in the interior of the global domain, and 
are 

M 

l N*^ - 2M 

j=l 

in number, where we recall that N** is the number of parameters in element . 

In the usual case where N* 1 is the same for all , one can specify exactly 
- 2. smoothness constraints at each of the (M - 1) interior nodes, leaving 
exactly N J - 2 conditions to be specified at the two ends of the global domain. 

If there are no additional end conditions on the function f(x) to be represented, 
the N** - 2 conditions must be arbitrarily specified and apportioned at the two 
ends. Such splines are therefore not unique. It is also clear that an even 
number for N* 1 will prevent a bias towards one end. While splines can be 
constructed for arbitrary <f>?(x), in most applications they are limited to poly- 
nomials. 

One can again construct cardinal basis functions 4>^(x), and employ represen- 
tation (41). Since the smoothness constraints couple the elements together, the 
cardinal basis functions are not at all localized, but extend over the global 
domain. They are thus inconvenient for computational purposes. There are two 
approaches that are used. In one, the original basis $?(x) is used, and the 
derivatives u . , u ., etc., are introduced as additional unknowns. The 
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interpolation and matching conditions enable one to solve for these derivatives 
* 

in terms of the u^, by inverting banded matrices. The other procedure, valid 
for equal intervals, is to introduce a new basis (x), known as B splines, 

which possess the smoothness property, but do not have the cardinal property 

B i 

(Xj) = <5^ . The B splines are non-zero only over N J elements, and thus 

have a localized nature. The coefficients of B spline expansions can again 

* 

be obtained in terms of the u^ by inverting banded matrices. A popular choice 
for polynomial splines is piecewise cubic (N- 1 = 4) , which leads to easily 
invertible tridiagonal matrices. An elementary discussion of splines is found 
in reference 1. 

Our discussion of functional approximations was aimed at their application 
in numerical solutions of operator equations. Their obvious role is to obtain 
expressions for derivatives and integrals in terms of nodal parameters, and 
to evaluate functions at points other than nodes. There are several other 
applications, which should be briefly mentioned. 

One application is to use piecewise functional approximations to obtain 
approximate analytic solutions of certain differential equations. In this 
method, known variable coefficients are replaced by simpler piecewise representa- 
tions in terms of known nodal parameters, so that the resulting equations 
possess an analytic solution in each element. The unknown solution coefficients 
are obtained by matching the solutions and their derivatives at interelement 
boundaries and applying boundary conditions at the global boundaries. The 
solution of the differential equations is thus reduced to that of an algebraic 
system for the coefficients. Further details are found in the works of Gordon 
(ref. 4) and Canosa and de Oliveira (ref. 5). 

Another important area of application is the representation of complex 
surfaces in physical space. The independent variables x are the two parameters 
defining parametric curves on the surface, while u stands for the position 
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vector. If the surface is very complex, one needs a piecewise representation, 
dividing the surface into patches. One class of such representation uses tensor 
products, interpolating through data given at the comers of the patches. 
Examples are programs developed at McDonnell Douglas (ref. 6), using piecewise 
cubic Hermite polynomials, and the work of Riesenfeld (ref. 7) employing B 
splines. In another class of representation, the curves defining the boundary 
of the patch are analytically prescribed, and one seeks what are referred to as 
blended interpolations for points inside the patch. Examples are the work of 
Coons (ref. 8) using Hermite polynomials, and Gordon (ref. 9) using splines. 

All of these approximate surface representations can play an important role in 
generating finite difference and finite element grids and formulating surface 
boundary conditions, for the solution of flows past complex boundaries. 

In closing we list the various degrees of freedom in a functional approx- 
imation. One is the choice of single versus piecewise representation, and 
the nature, size, and location of volume elements in the latter instance. 

Another is the functional form, which involves a choice of basis functions in 
the linear case. The number, location, and nature of interpolating nodes 
is another degree of freedom. Finally, for piecewise representations, there 
is the choice of using additional smoothness constraints to define splines. 

While continuity and convergence criteria make some of these choices inter- 
dependent, it still allows for large degree of flexibility in constructing 
functional approximations. 

FORMULATION OF THE EXACT EQUATIONS 

There are two mathematically equivalent ways to formulate the equations 
describing continuous fields. In the direct approach, the equations and 
boundary conditions of the problem are given. For certain classes of 
equations, an indirect variational formulation is possible, which incorporates 
some of the boundary conditions. A finite difference numberical solution is 
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usually based on the direct formulation, while the variational formulation is 
the starting point for one form of the finite element method. These two 
formulations are briefly discussed below, where x will stand for the complete 
set of independent variables. 

Direct Formulation 

The normal way to formulate a field problem is to specify that u(x) is a 
solution of an operator equation (s) 

G(u) = 0 . (58) 

A complete specification requires subsidiary equations valid on subspaces of 
x called boundaries. If S 1 represents the ith boundary subspace, the 
subsidiary equations 

B 1 (u) = 0 , xeS 1 (59) 

are called the boundary conditions on S 1 . The boundary S 1 can be prescribed 

or free, i.e., implicitly defined in terms of another operator equation. In 

many problems S 1 is defined in the limit as x approaches infinity. If 

S = Z S 1 defines a closed subspace, then (58) and (59) define a boundary value 
i 

problem. If x is one -dimensional, one can also have an initial value problem, 
where all the boundary conditions are specified at only one boundary. For 
multidimensional x, a mixed type of initial boundary value problem is possible, 
which is an initial value problem with respect to one (time- like) independent 
variable, and a boundary value problem in the subspace defined by the other 
independent variables. 

Two other points should be made with respect to a problem formulation. 

In certain problems, internal boundaries (such as shocks or slip surfaces) can 
occur, where the solution is discontinuous. The conditions at these surfaces 
are not boundary conditions in the sense used here, since they are actually 
limiting forms of the field equations (58). The other point refers to certain 
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classes of boundary value problems, in which a well behaved solution only 
exists for specific values of certain parameters, called eigenvalues. In 
eigenvalue problems, the determination of the eigenvalues can be an important, 
if not the principal objective. 

Variational Formulation 

An indirect variational formulation exists for boundary value problems 
in which the operator G(u) (58) is self adjoint. Then G(u) is a Frechet 
derivative of another operator F (u) , i.e., 

G - F u . (60) 

The operator F (u) defines the integral functional I(u) given by (20), whose 
Frechet differential is given by (22) . A variational statement of the original 
problem states that 

<5I(u,6u) = 0 (61) 

for all variations 6u. This immediately implies (58), the original operator 
equation, which is then referred to as the Euler equation. But it also 
requires the boundary conditions (59) on each S 1 to be such that the boundary 
integral in (22) is equal to zero. If only the first term exists, the require- 
ment is that either 6u is zero, i.e., u is prescribed on the boundary, or 

7 is zero. Thus (59) would be limited to 
o 

B 1 (u) = u - <J> 1 (x) (62a) 

or B 1 (u) = F q (u) , (62b) 

where cf> x (x) is a prescribed function of x on the boundary S 1 . If the second 
term also exists in the boundary integral in (22), we additionally require that 
either u n is prescribed, or 7^ is zero, etc. The conditions u, u n , etc., 
prescribed are called the principal boundary conditions, while the alternate 
conditions 7 q = 0, = 0, etc., are called the natural boundary conditions. 
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Thus, corresponding to each term in the boundary integral in (22), a boundary 
condition (59) must exist on each S X , which is either the principal condition, 
or the natural condition determined implicitly by (58). The variational 
statement (61) is thus subject to the constraints of the principal conditions, 
but automatically incorporates the natural conditions. 

There are problems for which G(u) is self adjoint, which involve boundary 
conditions (59) that are neither principal nor natural, as defined above. It 
is usually possible to extend the variational principle to include those 
cases. Let 

H 1 (u) = h 1 (x,u,u n ,u m , ...) , xes 1 (63) 

be a local operator defined on the boundary subspace S 1 . The Frechet 
differential can be written as 

H^u.du) = H*6u + H^<Su n + H^6u + ... , xeS 1 (64) 


Then the extended integral functional 


I(u) 


= J F (u)dx + l j H 1 


(u)dx 


(65) 


has the Frechet differential 

61 (u, u) = J F u 6udx + Z J [(F o + H*)6u + (Fj+H^fii^ + (F^H^du^ + ...]dx . 

V S 1 (66) 


The extended variational principle (61) now possesses extended natural boundary 

conditions ¥ + H 1 = 0, F, + H x = 0 , etc. For most cases, one can find 
o o ' l 1 

operators H 1 (u) such that boundary conditions (59) that are not principal 
conditions can be made to be extended natural conditions as defined by the 
extended functional (65). One can also show that several different choices 
for H 1 ^) are possible in some situations. 

When (58) or (59) involve several equations it is possible to handle some 
of them using Lagrange multipliers. Specifically, if (58) or (59) are 
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replaced by 


and 


G (u) =0 

and G q (u) = 0 


(67 

B 1 (u) = 0 

and B^(u) = 0 
0 

xeS" 1 , 

(68) 


where G(u) is the Frechet derivative of F(u), the variational principle can 
then be stated in terms of the functional 


I(u) 


■l 


[F(u) + XG q (u) ] dx + E 


[H 1 (u) + p 1 B^(u)]dx , 


(69) 


where the functions \(x) and y 1 (x) are parameters to be varied independently. 
G(u) = 0 and G q (u) = 0 are the Euler equations corresponding to the variations 
of <5u and 6A, respectively. Similarly, some of the equations B^u) = 0, and 
B q(u) = 0, are the natural conditions corresponding to the variation of 6u, 
6u n , ..., and <SA on the boundary S 1 . Sometimes the roles of G(u) and G q (u) can 
be reversed, leading to alternate variational principles for the same problem. 

The variational formulations discussed so far have been restricted to 
prescribed boundaries S 1 . We indicate the modification due to a free boundary 
by considering the case where the boundary S 1 and boundary condition B^u) 
are determined by the solution of 


g 1 [u(x) > x] = 0 , 


(70) 


which implicitly defines S^u) and V(u). The functional I (u) is now written 
as 


I(u) = 


J F(u)dx , 
V (u) 


(71) 


and its Frechet differential is 


6l(u,6u) = J <5F (u, 6u)dx + E | F(u)6n 1 dx 

V V 


(72) 
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where 611 1 is the amount "S’ 1 moves normal to itself due to the variation 6u(x), 
and V and S’ 1 are the domain and boundary before u is varied. If SnfS* 1 ) 
represents the variation 6u at the fixed boundary S’ 1 , one can show from (70) 


that 


.i 

6n = — = 


g + g u 
*n *u n 


(73) 


Combining (72) and (73), we find that the free boundary modifies the natural 
boundary condition to read 

pg, 3 : 

= 0 . (74) 


F “T 

o 11 

s + s u 
6 n 6 u n 


The variational formulation has two advantages over the direct formulation. 
The operator F(u) is a lower order operator than G(u), permitting a functional 
approximation with a lower degree of continuity. Also, since the variational 
formulation has the natural boundary conditions built into it, it therefore 
has fewer boundary conditions to satisfy than the direct formulation. It has 
the disadvantage of being indirect, and only existing for a certain class of 
problems. 

We are now ready to examine how the two types of discretizations discussed 
in the previous section are used to obtain approximate solutions to continuous 
field problems, starting from either of the two formulations discussed above. 


FINITE DIFFERENCE METHODS 

Any approximate method of solving a continuous field problem whose starting 
point is the discretization of some of the independent variables will be termed 
a finite difference method. The most common procedure employs point discretiza- 
tion at N nodes x^, with their associated unknown nodal parameters which can be 
functions of the variables t. This lends itself naturally to the direct 
formulation, by evaluating (58) approximately at N evaluation points Xy which 
do not necessarily coincide with the x^. (Recall that if u stands for several 
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dependent variables, each may be discretized by a different set of nodes.) The 
operator G(u) involves differential operators in x which must be approximated 
by finite difference operators in terms of the unknown nodal parameters. This 
has two consequences. In order to obtain simple difference approximations, 
it is highly desirable to choose the nodes to lie along coordinate surfaces if 
x is multidimensional. The other point refers to the nature of the functional 
approximation which is implied by the difference approximation. In the 
previous section, functional approximations were defined over global regions. 

Yet in initial value and initial boundary value problems, the solution along 
the time- like coordinate is only known up to the point presently reached in the 
calculation. Even in boundary value problems, a difference approximation based 
on a global functional approximation would be overly complex, and result in the 
need to invert very dense matrices. For these reasons, traditional finite 
difference approximations are based on functional approximations that interpolate 
data at nodes x^ limited to a neighborhood of the evaluation point x\ . We 
discuss such local difference approximations first, and subsequently examine 
some recent difference approximations based on global functional approximations. 

Methods Based on Local Functional Approximations 
The first observation one should make is that local functional approxima- 
tions used at neighboring evaluation points are in general incompatible. This 
can be simply seen by considering second order Lagrange polynomial interpolations 
in one dimensions for - x^ Using symmetrically placed points (leading to 

central difference formulas), the local functional approximation at is a 

* * * _ 

parabola through the points u^_^, » and u ±+i > while that at x^ + ^ is a 

* * * 

parabola through u i , u i+1 , and u i+2 . These two approximations describe two 

different curves in their region of overlap between x i and x i+1 * Once the 

* 

approximate solution for the uj is obtained, it is not clear which of the 

* 

curves to use in order to interpolate for the values of u between nodes 
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(presumably a weighted average would give the best results). Actually, the 


question is rather academic, since the difference in the values given by the 


two approximations should be no greater than the errors in the approximations 
* 

u. . 

J 


By contrast, the piecewise functional approximations of the previous 
section, although localized in nature, are disjoint functions that butt together 
with no regions of overlap. They therefore give unambiguous values for any 
quantity (except derivatives at interelement boundaries of an order higher than 
that demanded by the smoothness of the approximation) . Yet it is this very 
ambiguity in the local functional approximation which gives a local finite 
difference approximation its flexibility and power. If G (u) is quasi-linear, 
the local value of u determines the nature of the operator, which in tern 
determines the optimum type of difference approximation. Thus the nature of 
the local approximation can be determined at each evaluation point by the 
local solution. This is the basis for upwind differencing and the type 
differencing of transonic flows. Even at the same evaluation point, different 
terms in the operator G(u) can be approximated separately. The nature of the 
approximation can be made to change during an iterative solution, or a marching 
solution with respect to another independent variable. ADI methods and 
splitting or factorization techniques are applications of this degree of 
flexibility. 

The local functional approximation also has to be modified for evaluation 
points x . s near or at global boundaries, in order to satisfy boundary conditions. 
This can be done most readily if the global boundary is a coordinate surface. 

For a more general boundary which does not conform to the coordinate system, 
the approximation can become quite involved, if one wants to maintain the same 
level of accuracy. For this reason, a nonconforming boundary should be 
avoided if possible. 
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Lagrange Representation 


The simplest types of finite difference formulas for derivative operators 

are based on Lagrange polynomial interpolation. This is the basis for the 

standard forward, backward and central difference formulas for partial 

derivatives of any order, and of any order of accuracy. Lagrange interpolation 

* * 

can also be used to express u^ in terms of neighboring u^ (j ^ i), assuming 
* 

that u^ is unknown. Such a device is used in some numerical algorithms. 

In solutions involving time-like coordinates, final values of derivatives 
are already known at points previously computed. In boundary value problems, 
one needs to compute the same derivative at all nodal points. This suggests 
the use of Hermite interpolation to provide more accurate difference formulas 
without increasing storage requirements. 


Hermite Representation 

An example of a Hermite finite difference formula in one dimension is 

* * * * * 

derived from the specification of u^^, u xx (i-l) ,u i » u i+l and u xx(i+l )* whic ^ 

define a unique quartic polynomial. (This is an example of a Hermite representa- 

* 

tion with nonconsecutive derivatives.) Evaluating the second derivative of u (x) 
at = x^ (for equal spatial intervals h) , one obtains 

h 2 (u* n + 10 u* . + u* .J = 12 (u? + 2u* + u* ) , (75) 

V xx(l-l) XXI XX(l+l)' v 1-1 1 1+1 

which is the standard Hermite centered finite difference formula (ref. 10). The 

</c ^ 

solution for u xx ^ is obtained by tridiagonal inversion. Other Hermite differ- 
ence formulas involving any partial derivatives can be similarly obtained. 

An important application of Hermite interpolation is the construction of 
difference formulas for initial value problems of the form 


G (u) = u - f(x,u) = 0 , 


(76) 


where x is one- dimensional. If the solution is known up to the point x^, the 

values of u. and u . for all j < i are available to construct a variety of 
I xj 
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A more 


local Hermite interpolates from which one can explicitly predict ll +1 . 

accurate but implicit difference formula is obtained by including u x ^ + 2) 

the representation. Such a formula is normally used as a corrector in an 

iterative solution, where a predictor formula and (76) were first used to 

* 

calculate a first approximation for u x ^ + 2 )* 

Another approach to the numerical solution of initial value problems 
employs higher derivatives u , u , etc., which can be obtained in terms of 

Xa XXX 

lower derivatives by differentiating (76). One can then construct the Taylor 


* * * 1 * 2 
u (x) = U. + u . (x - X. ) + T U . (x - X. ) + . . . 

u ' l xi i J 2 xxi^ x J 


If the series is truncated after a finite number of terms, the result can be 
looked at as a local Hermite interpolate through the single point x^. Thus any 
step in a finite difference algorithm for the solution of an initial value 
problem can be obtained from a local Hermite interpolation (although the local 
functional approximation corresponding to a given algorithm is not necessarily 
unique) . 

Different representations can be used in obtaining difference formulas for 

initial boundary value problems. Let t be the time- like variable, and assume 

that by differentiating (58) one can express u^., u^, etc., as functions of 

u x* u xx* etc# > where x represents the remaining independent variables. If 

* 

one first discretizes x space, and defines Lagrange parameters u^(t) at the 

nodes x^, one can then use a local functional approximation and Lagrange inter- 

* 2*2 

polation to evaluate u x , u xx , etc., and obtain expressions for du^/dt, d u^/dt , 
etc. The latter can then be used to define a Hermite discretization of the t 
coordinate, and the solution can be advanced in t, using (77) (with x replaced 
by t), which represents Hermite interpolation through a single point in t space. 

In summary, any standard finite difference algorithm for solving a set of 
partial differential operator equations can be derived by applying sequences of 
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local functional approximations, and interpolating parameters of Lagrange or 
Hermite representations. The number of points and the order of interpolation 
determine the accuracy of the approximation (i.e., truncation errors). This 
still leaves freedom in the choice of points and parameters, and nature of the 
functional approximations. These can all be optimized to provide the best 
stability properties for the numerical solution. 

Methods Based on Global Functional Approximations 
We turn now to finite difference methods based on global functional 

approximations, limiting ourselves to Lagrange discretization at nodes x, and 

* 

the case x^ = x^ Thus the nodal parameters are u.^(t), where t represents the 

remaining undiscretized variables. Partial derivatives are special examples of 

linear operators obeying the property 

L(au + bv) = aL(u) + bL(v) , (78) 

where u and v are two arbitrary functions, and a and b are constants. Thus a 
local operator G(u) can be written generally as 

G(u) = g[x,t,u, L t u, L x u, L x (L t ijO , (79) 

where the subscripts indicate the variables on which the linear operator L 

operates, and g is an arbitrary function of the six arguments. For any set of 

linearly independent basis functions 4>^(x) the linear representation (28) can 

be expressed in terms of canonical basis functions 4 k (x) and the nodal 
* 

parameters u. (t) as 

1 * N * ^ 

u (x,t) Z u. (t) cf). Cx). (33) 

i=l 1 1 


The two basis functions are related by defining matrix elements a-j as 


a ij = W > 


(80) 


and the inverse matrix with elements b^ satisfying 


b ij a jk = 6 ik * 


(81) 
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Then 


(82) 


= b i;j <J>j(x) . 

If we define 

00 = L x [* ± CxD] , (83) 

it follows from (33), (82), and (83) that 

* N * 

L lu (x,t)] « Z u. (t?)G (x) , (84) 

x i=l 1 1 


where 

X^(x) = b ijXj(x) = L x [^.(x)] . 

For linear operators operating on t we obtain 

N 

L [u*(x,t)] = Z L (u*(t)]0.(x) . 
Z i=l x 1 1 


(85) 


( 86 ) 


Evaluating (79) at the evaluation points x^ = 

* 

of equations for the parameters u^(t): 

N 

g[Xj,t,u^ (t), L t [u*(t)], _Z utctOx^x..), 


X. 


3 


we obtain the following sets 


Z L [u*(t)])T i (x )] = 0. (87) 

i=l J 


For an arbitrary set of <Jk (x) , the matrices a^. are dense, and their inversion 
is inefficient. The practical use of (87) requires restrictions on the 
functions 4>^(x). Three such choices will be described, each leading to a 
practical finite difference method in x space. An arbitrary, independent 
method can be used in each case to perform the numerical solution in the t 
space. 


Finite Fourier Series 

If x is one- dimensional, with periodic boundary conditions, a convenient 
choice is 

<|> k (x) = e 2irikx/L ^ ( 88 ) 

where L is the length of the region, and i = /T The representation is the 
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finite Fourier series 


u*(x,t) = I C k (t0e 2lrikx/L , (89) 

k=-K 

where N = 2K+1. If x . are chosen to be equally spaced, the transformation 
* 

between u^ (t) and the c^(t) (corresponding to matrix multiplication by a^ and 
b^) is accomplished efficiently by using fast Fourier transforms (ref. 11). 

For linear differential operators, the functions x^GO as defined by (83) are 
just proportional to the 4^00, leading to further simplifications. Finite 
difference methods using (89) are referred to as pseudospectral (ref. 12) or 
"accurate space derivative" (ref. 13) methods. They can be extended to higher 
dimensions using tensor products. 

A Fourier series can be reformulated as an expansion in Chebycheff 
polynomials, based on the identity 

n k 

cos nx = £ a v cos x . 

k=0 

The nodes x^ are no longer equally spaced in the new x domain, but are located 
at the roots of the Nth order Chebycheff polynomial. Thus the basis for the 
accuracy of such a difference scheme is the same one that underlies Gaussian 
quadrature . 

Differential Quadrature 

The ideas behind the polynomial formulation of a Fourier method can be 
generalized to any set of orthogonal polynomials, with the nodes x^ again 
chosen at the roots of the Nth order polynomial. The matrix elements X^(*j) 
are easily calculated, using the properties of the orthogonal polynomials. The 
method, known as differential quadrature, is described in reference 14. 

Spline Differencing 

A third approach using global functional approximations is to use piecewise 
approximations, with nodes and evaluation points located on interelement 
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boundaries. In one dimension, a spline approximation is necessary, with 
smoothness constraints determined by the highest derivative present in the 
operator L x . As indicated previously, cardinal basis functions are not 
practical, and the original basis functions are employed, with the derivatives 
at the interelement boundaries as additional unknowns. These derivatives are 
related to the nodal parameters through banded matrix relations, rather than 
explicitly as in (84). For a cubic polynomial spline, the relation for 
first derivatives is identical to the one given by Hermite differencing. The 
second derivative relation differs from the Hermite formula, with (75) replaced 
by 

h^(2u 8u . + 2u ,0 = 12(u. , + 2u. + u. .) . (90) 

Sequences of one-dimensional differencing using ADI methods and splitting, are 
used in multidimensional problems. Further details on the use of spline 
differencing in the numerical solutions of partial differential equations are 
found in reference 15. 


Finite Volume Differencing 

Finite difference equations based on volume discretizations are often 
employed when the operator equation (58) can be written in divergence (or 
conservative) form 

& *4 


Here F(u) is a locator operator on u. Let x be discretized into N volume 
elements V 1 , each of which is enclosed by a set of boundaries S 1 ^. If (91) 
is integrated over element V 1 , and the divergence theorem (4) is applied, the 
result can be written as 


auM + \ z s ij F ij = o , 

V 1 j 


3t 


where u** is defined by (25) , and F ^ is defined as 


(92) 
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★ 

FCu ) ndx . 


(93) 


The unknown parameters are the u , and a functional approximation is required 
to express the average normal fluxes F 1 ** in terms of these parameters. This 
can be most easily demonstrated for one-dimensional differencing. Letting the 
subscript k refer to a global numbering of interelement boundaries (which are 
nodes in one dimension) , a local functional approximation for element V 1 can 
be written as 

u 1 (x) !=> I U k ^(x) . (94) 

k (i) 


The notation k(i) indicates a particular choice of nodes k in the neighborhood 
of element V 1 , and^(x) is a local canonical basis function. The latter is not 
to be confused with the element cardinal basis functions defined by (43) and 

i ^yi 

(44). The nodes k(i) need not be contained in V , and ^ (x) ^ 0 for those 
nodes. Integrating (94) over element V 1 , one obtains 

u* 1 = E U*^ • (95) 

k(i) K k 


where 

l5 Cx)dx • 

V 1 


(96) 


* *■[ 

The u^ can then be expressed in terms of the u by inverting a sparse matrix 

in a manner similar to that which exists for Hermite differencing. 

The determination of the F 1 * 1 , when F(u) involves differential operators, 

again creates ambiguities resulting from the incompatibility of local 

i * 

functional approximations at neighboring elements V . Once the are 

* i 

obtained by inverting (95) , the local representation u (x) can be obtained 

*i 

from (94) . One can then determine F (u ) , and use (93) to calculate the 
terms F 1J for the two boundaries along the x direction. If this procedure is 
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followed, the value of F 1J for a given boundary separating two elements will 
be independently calculated for each of the elements. Yet global conservation, 
obtained by summing (92) over all the elements, requires that the two F 1 * 1 be 
equal in magnitude and opposite in sign. One must therefore choose a single, 
unambiguous F 1J for each boundary, using some averaging or biasing. If t is 
a time-like variable, the bias can be alternated with each marching step. 

Note that at global boundaries exact prescribed values of F 1 ^ can be imposed. 

If x is multidimensional, the above one-dimensional differencing can be 
used sequentially along several coordinates, using splitting techniques. A 
particular advantage of finite volume differencing is that the original 
equation (58) has been integrated, so that the operator F (u) involves lower 
order differential operators than G(u). Therefore a cruder local approximation 
(94) can be employed. The possibility of alternating the bias when t is time- 
like, allows even still cruder approximation for each marching step (ref. 16). 
The conservative, or integral nature of the numerical solution also guarantees 
that jump conditions across discontinuities are automatically satisfied, even 
if the discontinuities are smeared out by the calculation. 


Methods Based on a Variational Formulation 

We conclude this section by describing briefly a finite difference 

approach based on the variational formulation (61) and (65). The starting 

point is the same as for the direct formulation. One first chooses a set 

of nodes x.. and evaluation points x. inside the domain V and on the boundaries 

S 1 , the type of nodal representation, and the nature of the functional approxi- 

* i * 

mation. These are then used to evaluate the operators F (u ) and H (u ) at 
the evaluation points x^ , as functions of the nodal parameters. The next step 
is to approximate the integral functional I (u) in terms of the discretized 
F (u) and H 1 (u) . This is done by appropriate quadrature formulas of the form 
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( 97 ) 


I(u*)^Z w, [F (u*) ] Zwj [H^u*)]. 

j 3 3 i k K K 

Here w^ and w£ are weight coefficients defined implicitly by some functional 

approximations of the respective integrands, (These functional approximations 

* 

can in general be independent of those used in obtaining [F(u )]. and 
i * 

[H (u )]j in terms of the nodal parameters.) The summations in j and k are 

over the evaluation points contained in the domain V and boundary S 1 y 

respectively. In many cases one simply chooses w. = w? = 1. 

J K 

* * 

With I(u ) expressed as a function of the nodal parameters u^ through 

(97) (assuming Lagrange representation for the moment), the variational 

principle (61) simply becomes 

— — y = 0 , j = 1 to N , (98) 

3u. 

J 

providing N equations for the N unknown parameters. This method is sometimes 
called the Euler method, and is described further in reference 17. Actually, 
the method bears a striking resemblence to methods based on functional 
approximations, being somewhat hybrid in nature, with one foot in each camp. 

It is therefore a good point to leave finite difference methods, and turn out 
attention to functional approximation methods. 

FUNCTIONAL APPROXIMATION METHODS 

Any approximate method of solving a continuous field problem whose 
starting point is the discretization of the dependent variables will be 
termed a functional approximation method. We will describe such methods in 
terms of the general functional representation (26), applying the approxima- 
tion first to a variational formulation, and subsequently to the direct 
formulation. The results will then be specialized to linear and piecewise 
representations, the latter giving what we normally called finite element 
methods. 
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Variational Formulation 

The functional approximation (26) lends itself naturally to the varia- 
tional fonnulation (61) and (65). The method will be first described for the 
case when the dependence on all the variables will be approximated so that 
the functions c^ in (26) become constant, and there is no variable t. This 
case is usually called the Ritz or Rayleigh-Ritz method. The function 
g(x;Cj) must first be chosen so as to satisfy the principal boundary 
conditions. Substitution of (26) into (65) yields 


I(Cj) ^ j F[g(x;Cj)]dx + S 


H 1 [gCx;c j :)]dx . 


(99) 


This restricts g(x;c^) further to functions with sufficient continuity for 
the integrals to exist. The variational principle (61), applied to all varia- 
tions $Cy gives the set of equations 

ei 0 j = 1 to N (100) 


for the N parameters c^ . 

The method can be extended to functional approximations (26), where c^. 
are now functions of undiscretized variables t. It is then referred to as 
the Kantorovich method. The integral functional (65) must now be written as 


I(u) 


j ^ F (u) dxdt + Z | | H 1 (u)dxdt , 

st aht) 


( 101 ) 


T V(t) 


where T and refer to the subdomain of variables t, and their boundaries. 
Substitution of (26) into (101) yields 


Uc. (t)] 

T 


J F{g[x-c.(t)]}dx]dt + Z J [ 
V(t) sj 


^ H^gtxjc.. (t)]}dx]dt. 

s A (t) 


( 102 ) 
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Equation (102) is now considered an integral functional over t space involving 
the unknown functions Cj(t). The variational principle (61) then states that 
the Frechet differential 

fiKcj.acj) = o , ( 103 ) 

which results in the set of equations and boundary conditions necessary to 

determine the set of unknown functions c.(t). 

J 

As indicated before , the operator F (u) involves lower order differential 
operators than G(u), permitting a functional approximation with lower order of 
smoothness. The approximation also need not satisfy the natural boundary 
conditions, since they are automatically satisfied in the variational process 
(to the same degree that the equation G (u) = 0 is satisfied). For these 
reasons a Ritz or Kantorovich method is much to be preferred. Unfortunately, 
it is limited to boundary value problems in which G(u) is self-adjoint. There 
have therefore been many attempts to create so-called "variational" principles 
designed to solve problems for which a true variational principle does not 
exist. These new principles may be classed as adjoint variational, quasi- 
variational, or restricted variational. Finlayson and Scriven (ref. 18) have 
shown that they are all either based on a direct formulation in disguise, or 
offer no real advantage over a method based on a direct formulation. There 
is therefore no further need to consider any of these formulations. 

A new method which makes use of a variational formulation in an iterative 
procedure is the pseudo- functional method of Norrie and deVries (ref. 19). 

It is designed for problems which come close to admitting a variational 
principle. More precisely, assume that (58) is given by 

G(u) = F u (u) + G o (u) = 0 , (104) 

and the boundary conditions (59) that are not principal conditions can be 
written as 
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( 105 ) 


B X (u) = F.(u) + B*(u) = 0 , xeS 1 , j = 0,1, etc., 

where the operators F\ (u) are related to F(u) as in (20) and (22). If the 
terms G q (u) and B^(u) are sufficiently small, then (104) and (105) can be 
solved by iteration. Let u m * represent the solution after m-1 iterations. 
Then u m is defined as the solution of 

F u (u* m ) + G o Cu* (m-1) ) ^0 , C 


subject to the natural boundary conditions 

T. (u m ) + B ^ (u Cm ’ i;) ) ^0, xeS 1 , j = 0,1, etc., (107 

Equations (106) and (107) are thus seen to follow from the application of the 
variational principal (61) to the functional 

I(u* m )^ [ [F (u* m ) + G (u* (m ' 1) )u* m ]dx + £ f H 1 (u* (m " 1) )u* m dx . (108 


An iterative Ritz procedure can be applied to (108), until a converged solution 
for the Cj is obtained. 

The Method of Weighted Residuals 

If a variational formulation does not exist, even approximately, then a 

functional approximation method must be based on a direct formulation. To 

accomplish this, (58) and (59) must be converted into functionals. To see 

how this can be done, let us rewrite the Ritz procedure applied to a variational 

formulation, in terms of the equivalent direct formulation. If we substitute 

(26) and (66), and apply the variational principle (61) to all variations 

6c. , we obtain 
J 


fc? G < u *> dx ♦ l jiff? B >*> + 3?r 


+ 7 ^ nn B^ (u ) + ... ]dx ** 0 , (109) 

j 
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where we used (60) and let = F q + = 0, B* = + H* = 0, etc., repre- 

sent the extended natural boundary conditions in (59) . (The principal 
boundary conditions give 3g/8c.. = 0, 8g n /8c.. = 0, etc., on S 1 , and are assumed 
satisfied by the choice of g(x;c^). Thus contributions to the boundary 
integral terms in (109) will only come from the natural boundary conditions.) 

By using integration by parts and the divergence theorem (4), one can show the 
equivalence of the sets of equations (109) and (100). But (109) could have 
been obtained from the direct formulation by integrating (58) and the natural 
boundary conditions (59) over their respective domains, after first multiplying 
by appropriate weighting functions. Particular linear combinations of these 
integrals then yield (109). Note that the weighting function for B^(u) is the 
same as that for G(u), but those for B*(u), B^u), etc., (if they are present) 
are different. 

The above considerations suggest that the direct formulation (58) and (59) 
be recast in the equivalent weak form 

J tKx) G(u)dx = 0 , (110) 

V 

and 

<Kx) B 1 (u)dx = 0, (111) 

S 1 

where (110) and (111) are assumed valid for all arbitrary functions i|>(x). A 

functional approximation method can be obtained by choosing a finite set of 

linearly independent weighting functions ^ (x) to approximate ^(x), and sub- 

* 

stituting (26) and each ^ (x) in turn into (110). If G(u ) is termed the re- 
sidual, the resulting set of equations is thus obtained by equating to zero the 
integrals of weighted residuals over the domain. The method is therefore often 
referred to as the method of weighted residuals. If all the boundary conditions 
are not satisfied by the choice of (26), additional boundary residual equations 
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(Ill) must be calculated. These normally use the same weighting functions as 
in (110), although (109) shows that different weighting functions may be 
appropriate for some B X (u). 

By analogy with the variational case, integration by parts can be used to 
obtain integrals involving lower order differential operators. It is also 
possible to combine equation residuals and boundary residuals in the same 
equation, as was done in (109). To indicate these procedures, consider a 
term in G(u) that can be written as a divergence 8F/8x. Then the integral 
for that term can be written as 

j *3 H dx = j ^j nFdx " f ^jx Fdx * < 112 ) 

V S V 

If one of the terms in B x (u) is nF(u), it is then clear how (110) and (111) can 
be combined to eliminate that term. Note that (112) imposes smoothness condi- 
tions on (x) . We will henceforth examine the method of weighted residuals 
based on (110) with the understanding that these can be transformed by 
integration by parts and combined with (111) to eliminate certain boundary 
residual terms. When this is not possible, boundary integrals (111) would be 
treated in the same manner as (110). 

Let us generalize (110) by introducing the undiscretized variables t, 
and considering integrations over the domain and boundary of the discretized 
variables x. Thus, given (58), (26) and a set of weighting functions if^(x,t), 
the method of weighted residuals gives the equations 

I Gig (x;c^ (t) ]dx 0 j » 1 to N (113) 

V 

for the unknown functions c^ (t) . There are many possible choices for ij^(x,t), 
each one leading to a different method. They are fully discussed in the book 
by Finlayson (ref. 20). The various classifications are briefly summarized below. 
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Method of Moments 


If ipj (x,t) form an arbitrary, linearly independent set of functions, we 
have the general method of moments. Normally, it is restricted to functions 
of x only, and are typically members of a complete set of functions. A 
popular choice is polynomials in x. 

Galerkin Method 

In the Galerkin method, the weighting function is chosen to give the 
same equations as those provided by a variational formula. It follows from 
(109) that we must have 

V x,t) = % > (114 > 

where g is considered a function of x and c^ in performing the partial deriva- 
tive, i.e., t is considered a fixed parameter. This is probably the most 
popular method, particularly in finite element applications. 

Least Squares Method 

In this method we set 

V x>t) = Hr teC x s c iCt)l » (115) 

where again G(g) is considered a function of x and c^. The name of the method 
becomes obvious on substituting (115) into (113) and interchanging integration 
and differentiation, to obtain 

S §7 | G 2 r g (x; Ci Ct))]dx^0 . (116) 

While (116) minimizes the integrated square of the residual, a more logical 

2 

procedure would be to determine the maximum value of G in the domain for a 
given choice of c^, and to minimize this maximum among all choices of c^. 

While this has been used historically, it is difficult to apply in practice, 
and has been superseded by (116). One disadvantage of the least squares method 
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L 


is that the order of the differential operators cannot be lowered through 
integration by parts. 


Collocation Method 

If we admit discontinuous functions for ^ , several new methods are 
available. Let (j = 1 to N) be a set of n arbitrary points in the x domain, 
called collocation points. Then if 


M = 6(x - Xj ) , (117) 

where <5 represents the Dirac delta function, substitution of (108) in (103) 
yields 

js« 0 , (118) 

i.e., the residual is set equal to zero at the collocation points; hence, the 
name collocation method. Note that integration by parts is not possible. 


Subdomain Method 

If one divides the domain V into arbitrary subdomains V , one can define 
a less violent alternative to the Dirac delta function; namely, the 
characteristic function 


ip. (x) = 1 
= 0 


if xeV^ 
if x*V j 


(119) 


Equation (113) now becomes 

| G[g(x;c i (t)]dx * 0 . (120) 

V j 

Thus the integrated residual is set equal to zero in each subdomain V* 1 ; hence, 
the name subdomain method. Note that terms in G(u) that can be written as a 
divergence can be converted via the divergence theorem (4) to integrals over 
the boundaries of V ^ involving lower order operators. This method is some- 
times called the method of integral relations. 
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Least Squares Collocation Method 

The methods described above can often be combined. An example is the 
least squares-collocation method. We start with the least squares method 
(106), and approximate the integral by appropriate quadrature formulas over 
M arbitrary collocation points where M > N. The resulting equations for 
the parameters c^ are 


a M , M 

w i“ [g( x i > c k (t)0] *2 w.GCgCx.j^Ct))] 


( 121 ) 


In practice, one often chooses w^ = 1. As M approaches infinity, the method 
approaches the least squares method. On the other hand, if M = N and 
3G[g(x. ;c, (t)]/3c. is non-singular, then (111) is reduced to the collocation 

IK j 

method (108) (assuming that all w^ are non-zero). 

The equivalence of the N-point quadrature approximation to the least 
squares method and the collocation method can be generalized to any residual 
method involving continuous weighting functions. Omitting the dependence on 
t, we can write the N-point quadrature approximation to (113) as 
f N 

J ^(x)G[g(x;c k )]dxs* w.^ fx i )G[g(x i 5c k )] 0 . 


( 122 ) 


This reduces to the collocation method if ^ (x^) is non-singular and w^ are 
non- zero. Thus, if the integrals in a residual method are too complex to 
evaluate analytically, and no integration by parts is employed, an N-point 
quadrature approximation is identical to a collocation method. By a judicious 
choice of collocation points x^, this method gives results whose accuracy is 
consistent with the original functional approximation. The choice can be 
made rational if the functional representation is linear, which is the case 
we consider next. 
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Linear Functional Approximation 

All the methods described so far have been considered for a general 
functional approximation (26). In practice, one normally uses the linear 
representation (28). This simplifies some of the methods. For example, the 
weighting function for the Galerkin method becomes 

^ (x) = 4 to , (123) 

i.e., the weighting functions are the basis functions themselves. If (x) 
are given by (78), i.e., if (28) is a finite Fourier series, then the Galerkin 
method’ is called a spectral method (ref. 21). 

A linear representation allows the introduction of nodal parameters as 
unknowns by choosing an arbitrary point discretization x^. It is then 
possible to establish a correspondence with finite difference methods. The 
most obvious one is through the collocation method. If the collocation points 
are identified as the evaluation points, it is evident that the collocation 
method is identical to a nodal finite difference method employing a global 
functional approximation. The method of differential quadrature has its 
analogue in the collocation method, where it is referred to as orthogonal 
collocation. 

Most conventional finite difference methods employ local functional 
approximations. Since those functional approximation methods based on discon- 
tinuous weighting functions (i.e., collocation or subdomain) yield equations 
evaluated at disjoint points or subdomains, one can generalize them by per- 
mitting local functional approximations. One can then say that all convention- 
al finite difference methods are collocation methods using local functional 
approximations. Similarly, one can consider finite volume difference methods 
as subdomain methods (with the divergence theorem applied) using local 
functional approximations. Finally, the Euler difference method may be 
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thought of as a variational -col location method using local functional approxi- 
mations. 

Weighted residual methods using continuous weighting functions, and 
which do not employ quadratures, can only be formulated in terms of a global 
functional approximation. Even for discontinuous weighting functions, a 
global functional approximation may be preferred. For complex domains, the 
integrals resulting from such a global approximation could not be calculated 
analytically. Even for one-dimensional or tensor product approximations, 
global functional approximations would lead to dense matrices. Both of 
these difficulties can be avoided by using piecewise functional approximations, 
which will now be discussed. 

Finite Element Methods 

Any functional approximation method using a piecewise functional repre- 
sentation is termed a finite element method. Thus, the domain of x is divided 

k 

into M volume elements V , called finite elements, and a set of N global 

* 

nodes x . , and their associated nodal parameters u^ft). (We assume a Lagrange 

representation for now.) For each element V we have a set of local coordin- 

k k k *k 

ates x , N local nodes x^, the associated local nodal parameters u^ (t) , and 

k 

element cardinal basis functions (x ). The latter are called element shape 

* k 

functions. The representation of u (x,t) in element V is 

N k 

u K (x K ,t) ** I U* K (t)8?(x K ) . (49) 

JU1 z 

The use of (49) in the two types of functional approximation methods will be 
briefly outlined. 

Variational Finite Element Method 

We will describe the Ritz method for simplicity. The variational formula- 
tion (99) and (100) can be easily reformulated in terms of the finite elements. 

k k 

Let I be the contribution to the integral functional (99) from element V , 
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given by 


■XM p i X “’ k ^ c ’ tk)i 

v k 

(124) 

k k 

Here 3x/3x represents the transformation Jacobians for elements V and the 

ki 

element boundaries b . The boundary integrals exist only for elements lying 

on the global boundary, with contributions coming from those boundaries S 1 

k k *k 

that border element V . From (126) one can then determine 81 /3u^ . By 

means of the mapping between local and global node numbers, this can be 
k * 

rewritten as 8 1 /du. , in terms of global nodal parameters. The variational 
principal (100) is obtained by summing over all the elements, i.e.. 



91 ^ 

k=l 3u. 
J 


M 

Z 


0 , j = 1 to N 


(125) 


Note that contributions to (125) come only from elements containing global node 
Xj, and the resulting equation involves only the nodal parameters contained 
in those elements. This insures sparse matrices in the solution of the 
algebraic syfctem (125). 


Residual Finite Element Methods 

The method of moments does not provide a useful finite element method, 
since the weighting function is not localized. The most popular method is the 

ft 

Galerkin method. Omitting the dependence on t for the moment, if x^ is the 

k 

local node in element V corresponding to global node x^ , it follows from 

k 

(123) and (113) that the contribution to (113) from element V is 

f'?{c* k ) G [ ? U ;X(x k )] ^-dx k . ( 126 ) 

Am J l m=l m m 9x k 

V K 

k 

If (126) is renumbered with a global node numbering, and written as 1^ in 
terms of global nodal parameters, then (113) becomes 
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M k 

Z I t^O, j = 1 to N , (127) 

k=l J 

where contributions to (127) again come only from elements containing node 
x ^ . If part of G(u) had been integrated as in (112) to create boundary 
integrals, then additional boundary terms would be needed in (126) for nodes x 
on the global boundary. 

The least squares finite element method is formulated in a manner 
similar to the variational method, based on (116). The collocation finite 
element method follows directly by substitution into (118). If some of the 
collocation points x. lie on element boundaries, then Hermite or spline 
representations are required to provide sufficient smoothness to calculate 
the operator G. Lower order representations are sufficient if all the colloca 
tion points are in the interior of elements. 

An important advantage of finite element methods is that prescribed 
boundary conditions on global boundaries are simply satisfied by setting the 
appropriate nodal parameters equal to their boundary values. Equations (124) 
or (126) would not be calculated for those nodes. Derivative boundary condi- 
tions can be satisfied by using Hermite shape functions. The number of un- 
known nodal parameters can be further reduced when some of the elements 

contain interior nodes. If x^ is an interior node located inside element 

k * 

V , then (125) (or (127)) is the only equation involving u. . The set of 

k 

equations for all the interior nodes in V can be solved for the interior 
nodal parameters in terms of the nodal parameters on the element boundary. 

By this process, called condensation, the final set of equations contains 
only nodal parameters associated with nodes on interelement boundaries. 

CONCLUSION ' 

Finite difference methods have been discussed from a rather unorthodox 
viewpoint in order to bring out their relationship to functional approximation 
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methods. Let us now examine this relationship by first comparing the nodal 
finite difference method with the finite element method. Both methods rely 
on a discretization of the x domain into nodes and the introduction of 
associated nodal parameters to represent the unknown function u(x). It is 
in the manner in which one obtains equations to solve for these parameters 
that the two methods diverge. 

The finite element method requires two additional discretizations. One 
is the discretization inherent in the global functional approximation which 
permits an unambiguous evaluation of u, or any operator on u, at an arbitrary 
point x. The other discretization, which is peculiar to the finite element 
method, is the additional discretization of the x domain into volume elements 
that define a piecewise functional approximation. These three discretizations 
are not independent, but are interrelated to provide desired smoothness to the 
approximation with the minimum of complexity. It is the achievement of these 
two contradictory goals that is the hallmark of the art of the finite element 
method. Finally, a variational or weighted residual method must be chosen 
to define appropriate integral functionals. The latter choice also involves 
some ingenuity, since integration by parts for continuous weighting functions, 
or proper choice of collocation points, can lessen the smoothness requirements. 
The choice of method is thus also coupled to the three discretizations. 

The conventional nodal finite difference method is essentially a 
collocation method, with nodes and collocation points aligned along coordinate 
lines if x is multidimensional. The finite difference approximations to the 
governing equations can be interpreted as resulting from local functional 
approximations. The approximation can therefore vary from point to point, 
and even for individual terms in equations. The art in the finite difference 
method is to use this great flexibility to obtain efficient and stable 
solutions. 
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The power of the finite element method lies in its ability to handle 
complex boundaries through the freedom in choosing the volume discretizations, 
and the ease in satisfying boundary conditions. An additional advantage 
exists for variational and certain weighted residual methods, where we can 
deal with operators of lower differential order and admit approximations of 
lower smoothness. Since a single global functional approximation is required, 
the method appears to be less flexible in dealing with the complex physical 
phenomena associated with highly nonlinear equations, such as those of fluid 
dynamics. Some progress has recently been reported in simulating type 
differencing (ref. 22) and upwind differencing (ref. 25) within the finite 
element method. 

The nodal finite difference method has the flexibility to cope with the 
phenomena associated with the complexities of the equations. On the other 
hand, boundary conditions can be satisfied accurately in practice only if the 
boundaries are coordinate surfaces. Here the recent work of Thompson (ref. 24) 
is generating coordinate systems for arbitrary surfaces gives promise to free 
the finite difference method from its major disadvantage. The use of piece- 
wise approximations (a finite element concept!) to represent arbitary surfaces 
can also play an important role. 

For initial boundary value problems, finite volume differencing can be 
thought of as the subdomain method with local functional approximations. Since 
it also results in lower order differential operators, it can be said to 
possess the other advantage attributed to finite element methods. Actually, 
its ability to treat discontinuities easily gives it somewhat of an advantage. 

In conclusion, finite element methods are best designed to handle complex 
boundaries, while finite difference methods appear to be superior for complex 
equations. Time and further research will tell if one of the methods will be 
able to overcome its shortcomings and emerge as clearly superior in solving 
boundary and initial boundary value problems. 
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